/**
 * Implement the intersection for 2D primitives. 
 * @author Songgang Xu (sxu at tamu dot edu)
 * @file intersect_2d.h
 */

#ifndef _INTERSECT_2D_H_
#define _INTERSECT_2D_H_

#include "predicates.h"



/**
 * Return the orientation of three 2D points. 
 */
template <typename Value>
int orient_2d (const Vector<Value> & p1, const Vector<Value> & p2,  const Vector<Value> & p3)
{
	Value v =  ((p1[0]-p3[0])*(p2[1]-p3[1])-(p1[1]-p3[1])*(p2[0]-p3[0]));	
	if (v > Value(0)) return 1  ; 
	else if (v < Value(0)) return -1 ;
	else return 0; //if (v == Value(0)) return 0 ; 
}


/**
 * Assume the two 2D segments are parallel. Compute the interseciton. 
 */
template<typename Value>
bool intersect_seg_seg_2d_parallel (const Segment<Value> & seg1, const Segment<Value> & seg2, IntersectResult<Value> & _intr)
{
	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
//	const Vector<Value> & d = seg2.p1;

	if (0 != orient_2d (a,b,c)) return false; 

	Segment <Value> s1 = seg1 ;
	Segment <Value> s2 = seg2 ;

	s1.order (); s2.order(); 

	if (s1.p0 < s2.p0)
	{
		//s1p0, s1p1 , s2p0, s2p1
		//if (s1.p1 < s2.p0) {_intr.pos = NONINTERSECT; return false;}
		//s1p0, s1p1 -- s2p0, s2p1
		if (s1.p1 == s2.p0 ) { _intr.res_type = POINT; _intr.point = s1.p1 ; return true; }
		//s1p0, s2p0, s1p1 s2p1
		if (s1.p1 > s2.p0 && s1.p1 <= s2.p1) { _intr.res_type = SEGMENT; _intr.segment = Segment<Value>(s2.p0, s1.p1) ; return true; }
		//s1p0, s2p0, s2p1, s1p1
		if (s1.p1 > s2.p1) { _intr.res_type = SEGMENT; _intr.segment = s2 ; return true; }
	}
	else if (s2.p0 <= s1.p0)
	{
		//s2p0, s2p1 , s1p0, s1p1
		//if (s2.p1 < s1.p0) {_intr.pos = NONINTERSECT; return false;}
		//s2p0, s2p1 -- s1p0, s1p1
		if (s2.p1 == s1.p0 ) { _intr.res_type = POINT; _intr.point = s1.p0 ; return true; }
		//s2p0, s1p0, s2p1 s1p1
		if (s2.p1 > s1.p0 && s2.p1 <= s1.p1) { _intr.res_type = SEGMENT; _intr.segment = Segment<Value>(s1.p0, s2.p1) ; return true; }
		//s2p0, s1p0, s1p1, s2p1
		if (s2.p1 > s1.p1) { _intr.res_type = SEGMENT; _intr.segment = s1 ; return true; }
	}
	return false; 
}

/**
 * predicate the point in convex polygon or not. 
 */
template<typename Value>
bool predicate_pol_convex_point (const Polygon <Value> & _pol1, const Vector <Value> & _point)
{
	unsigned size = _pol1.size(); 
	for (unsigned i = 0; i < size; ++i)
	{
		if (orient_2d(_pol1[i], _pol1[ (i+1) % size], _point) > Value(0)) return false; 
	}
	return true; 
}

/**
 * Predicate the interseciton of 2D segments
 */
template<typename Value>
bool predicate_seg_seg_2d (const Segment<Value> & seg1, const Segment<Value> & seg2)
{
	Value d0 = orient_2d(seg1.p0, seg1.p1, seg2.p0);
	Value d1 = orient_2d(seg1.p0, seg1.p1, seg2.p1);

	Value d2 = orient_2d(seg2.p0, seg2.p1, seg1.p0);
	Value d3 = orient_2d(seg2.p0, seg2.p1, seg1.p1);

	return (d0 * d1 <=0 && d2 * d3 <=0);
}



/**
 * Compute exact intersection of two 2D segments. 
 */
template<typename Value>
bool intersect_seg_seg_2d (const Segment<Value> & seg1, const Segment<Value> & seg2, IntersectResult<Value> & _intr)
{
	Value  s, t;       /* The two parameters of the parametric eqns. */
	Value num, denom;  /* Numerator and denoninator of equations. */


	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
	const Vector<Value> & d = seg2.p1;

	bool retValue = false;
 
	const int X = 0; 
	const int Y = 1; 
	denom = a[X] * ( d[Y] - c[Y] ) +
		b[X] * ( c[Y] - d[Y] ) +
		d[X] * ( b[Y] - a[Y] ) +
		c[X] * ( a[Y] - b[Y] );

	/* If denom is zero, then segments are parallel: handle separately. */
	if (denom == Value(0.0))
	{
		return intersect_seg_seg_2d_parallel(seg1, seg2, _intr);
	}

	num =    a[X] * ( d[Y] - c[Y] ) +
		c[X] * ( a[Y] - d[Y] ) +
		d[X] * ( c[Y] - a[Y] );
	if ( (num == Value(0.0)) || (num == denom) )  retValue = true; 
	s = num / denom;


	num = -( a[X] * ( c[Y] - b[Y] ) +
			 b[X] * ( a[Y] - c[Y] ) +
			 c[X] * ( b[Y] - a[Y] ) );
	if ( (num == Value(0.0)) || (num == denom) )  retValue = true; 
	t = num / denom;


	if      ( (Value(0.0) < s) && (s < Value(1.0)) &&
			  (Value(0.0) < t) && (t < Value(1.0)) )
	{

		retValue = true; 
	}
	else if ( (Value(0.0) > s) || (s > Value(1.0)) ||
			  (Value(0.0) > t) || (t > Value(1.0)) )
	{
		retValue = false;
	}

	_intr.res_type = POINT ;
	_intr.point = Vector<Value> (a[X] + s * ( b[X] - a[X] ), a[Y] + s * ( b[Y] - a[Y] ));

	return retValue;
}

/**
 * Justify if a point is in convex polygon. 
 */
template <typename Value>
bool is_in_polgyon_2d_cov (const Polygon <Value> & _pol, const Vector <Value> & _p)
{
	// Map the polygon and point to the 2D domain. 
	unsigned size = _pol.size();
	const unsigned x = 0, y = 1;

	for (unsigned i = 0, j = size - 1; i < size; j = i ++)
	{
		if (orient_2d(_pol[j], _pol[i], _p) < 0 ) return false; 
	}
	return true; 
}


template <typename Value>
bool is_in_polgyon_2d (const Polygon <Value> & _pol, const Vector <Value> & _p)
{
	// Map the polygon and point to the 2D domain. 
	unsigned size = _pol.size();
	const unsigned x = 0, y = 1;

	bool rvalue = false; 
	for (unsigned i = 0, j = size - 1; i < size; j = i ++)
	{
		if (_p == _pol[i] ) return true; 
		if ((((_pol[i][y] <= _p[y]) && (_p[y] < _pol[j][y])) ||
			 ((_pol[j][y] <= _p[y]) && (_p[y] < _pol[i][y]))) &&
			(_p[x] < (_pol[j][x] - _pol[i][x]) * (_p[y] - _pol[i][y]) / ( _pol[j][y] -  _pol[i][y]) + _pol[i][x]))
			rvalue = !rvalue;  					
	}
	return rvalue; 
}




template<typename Value>
int Advance( int a, int *aa, int n, bool inside, const Vector<Value> & v, std::vector<Vector<Value> > & result)
{
	if ( inside ) result.push_back(v);
	(*aa)++;
	return  (a+1) % n;
}

template<typename Value>
bool intersect_conv_poly_poly_2d (const Polygon <Value> & P, const Polygon <Value> & Q, IntersectResult<Value> & _intr)
{

	int inflag  = 0;
	const int Pin = 1, Qin = 2, Unknown = 0;  /* {Pin, Qin, Unknown}: which inside */

	int     a, b;           /* indices on P and Q (resp.) */
	int     a1, b1;         /* a-1, b-1 (resp.) */

	Vector<Value> A, B;           /* directed edges on P and Q (resp.) */
	int     cross;          /* sign of z-component of A x B */
	int     bHA, aHB;       /* b in H(A); a in H(b). */
	Vector<Value> Origin (Value(0.0), Value(0.0)); /* (0,0) */
	Vector<Value> p;              /* double point of intersection */
	Vector<Value> q;              /* second point of intersection */
	int     aa, ba;         /* # advances on a & b indices (after 1st inter.) */
	bool    FirstPoint;     /* Is this the first point? (used to initialize).*/
	Vector<Value> p0(2);             /* The first point. */
	bool     code;           /* SegSegInt return code. */ 
	int m , n ;

	std::vector<Vector<Value> > result ;

	/* Initialize variables. */
	a = 0; b = 0; aa = 0; ba = 0;
	inflag = Unknown; FirstPoint = true;
		n = P.size(); m = Q.size(); 

	do {

		/* Computations of key variables. */

		a1 = (a + n - 1) % n;
		b1 = (b + m - 1) % m;

		A = P[a] - P[a1] ; // SubVec( P[a], P[a1], A );
		B = Q[b] - Q[b1] ; // SubVec( Q[b], Q[b1], B );

		cross = orient_2d( Origin, A,  B );
		aHB   = orient_2d( Q[b1], Q[b], P[a] );
		bHA   = orient_2d( P[a1], P[a], Q[b] );

		/* If A & B intersect, update inflag. */
		IntersectResult<Value> intrseg; 
		code = intersect_seg_seg_2d( Segment<Value>(P[a1], P[a]), Segment<Value>(Q[b1], Q[b]), intrseg);

		if ( code == true && intrseg.res_type == POINT ) 
		{
			if ( inflag == Unknown && FirstPoint ) {
				aa = ba = 0;
				FirstPoint = false;
				p0 = intrseg.point ;
				//printf("%8.2lf %8.2lf moveto\n", p0[X], p0[Y] );
				//result_pol.push_back(Vector<double>(p0[X], p0[Y]));
			}
			p = intrseg.point ;
			// Add the intersection point here. 
			//std::cout<<p<<std::endl;
			result.push_back(p);
			if (aHB > 0) inflag = Pin; 
			else if (bHA > 0) inflag = Qin; 
		}

		/*-----Advance rules-----*/

		/* Special case: A & B overlap and oppositely oriented. They share one edge.*/
		if ( ( code == true ) && (( A* B ) < 0) )
		{
			if (intrseg.res_type == SEGMENT)
			{
				//std::cout<<"Shared Edge "<<intrseg.segment<<std::endl;
				_intr.res_type = SEGMENT ;
				_intr.segment = intrseg.segment ;
				return true; 
			}
			/*
			_intr.res_type = POINT; 
			_intr.point = p ;
			return true; 
			*/

		}
		/* Special case: A & B parallel and separated. */
		if ( (cross == 0) && ( aHB < 0) && ( bHA < 0 ) )
		{
			return false;
		}
		/* Special case: A & B collinear. */
		else if ( (cross == 0) && ( aHB == 0) && ( bHA == 0 ) ) {
            /* Advance but do not output point. */
            if ( inflag == Pin )
				b = Advance( b, &ba, m, inflag == Qin, Q[b], result );
            else
				a = Advance( a, &aa, n, inflag == Pin, P[a], result );
		}

		/* Generic cases. */
		else if ( cross >= 0 ) {
			if ( bHA > 0)
			{
				a = Advance( a, &aa, n, inflag == Pin, P[a], result );
			}
			else
				b = Advance( b, &ba, m, inflag == Qin, Q[b], result );
		}
		else /* if ( cross < 0 ) */{
			if ( aHB > 0)
			{
				b = Advance( b, &ba, m, inflag == Qin, Q[b], result );
			}
			else
			{
				a = Advance( a, &aa, n, inflag == Pin, P[a], result );
			}
		}

		/* Quit when both adv. indices have cycled, or one has cycled twice. */
	} while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) );
	
	if ( !FirstPoint ) /* If at least one point output, close up. */
	{
		//printf("%8.2lf %8.2lf lineto\n", p0[X], p0[Y] );
	    //result_pol.push_back(Vector<double>(p0[X], p0[Y]));

	}
	
	/* Deal with special cases.  */
	if ( inflag == Unknown)
	{
		if (is_in_polgyon_2d(Q, P[0]))
		{
			_intr.polygon  = P ;
			_intr.res_type = POLYGON; 
			return true; 
		}
		else if (is_in_polgyon_2d(P, Q[0]))
		{
			_intr.polygon = Q ;
			_intr.res_type = POLYGON ;
			return true; 
		}
		else return false; 
	}

	// There is a bug. The last and the first could be the same. 
    // Filter the duplicated vertices. 
    Vector<Value> first_point = result[0] ;
	unsigned pos = 0; 
	for (unsigned i = 1; i < result.size() ; ++i)
	{
		if (result[i] == result[pos]) continue; 
		result[pos + 1] = result[i] ;  
		pos++ ;
	}

	result.erase (result.begin() + pos + 1, result.end());
	if (result[result.size() - 1] == result[0]) throw ("Please fix the filter duplicated vertices bug");
	
	_intr.res_type = POLYGON ;
	_intr.polygon = Polygon<Value> (result);

	
	return true ;
}


/**
 - True. Same direction. 
 - False. Opposite direction.
 */
template<typename Value>
bool hseg_hseg_2d_parallel_dir (const Segment<Value> & seg1, const Segment<Value> & seg2)
{
	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
	const Vector<Value> & d = seg2.p1;
	
	return ( (b-a) * (d-c) > 0 ) ;
}

/**
 - Return Pin or Qin
 - Ture PIN
 - FALSE QIN
 */
template<typename Value>
bool hseg_hseg_2d_parallel_pqin (const Segment<Value> & seg1, const Segment<Value> & seg2)
{
	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
	const Vector<Value> & d = seg2.p1;
	

	if (a < b) {if  ( b <=d ) return true; else return false;}
	else {if ( b>=c ) return true; else return false; }
}
/**
 * Assume the two 2D segments are parallel. Compute the interseciton. 
 */
template<typename Value>
bool intersect_hseg_hseg_2d_parallel (const Segment<Value> & seg1, const Segment<Value> & seg2, IntersectResult<Value> & _intr)
{
	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
	const Vector<Value> & d = seg2.p1;

	// parallel lines
	if (0 != orient_2d (a,b,c)) return false;

	Value dot = (b - a) * (d - c); 
	if (dot < Value(0))
	{   // [ab] - [dc]
		if ( a <= c)
		{
			int count = 0; 
			if ( a >d && a <= c ) 
			{
				_intr.segment.p0 = a ;
				_intr.point = a;
				count++ ;
			}
			if ( c >=a && c < b)
			{
				_intr.segment.p1 = c ;
				_intr.point = c;
				count++ ;
			}
			if (count == 0) return false; 
			if (count == 1) {_intr.res_type = POINT ; return true;}
			if (count == 2) {_intr.res_type = SEGMENT ; return true;}			
		}
		// [cd] - [ba]
		else
		{
			int count = 0; 
			if (a >= c && a < d)
			{
				_intr.segment.p0 = a ;
				_intr.point = a;
				count++ ;
			}
			if (c > b && c <=a)
			{
				_intr.segment.p1 = c ;
				_intr.point = c;
				count++ ;
			}
			if (count == 0) return false; 
			if (count == 1) {_intr.res_type = POINT ; return true;}
			if (count == 2) {_intr.res_type = SEGMENT ; return true;}			
		}
	}
	else
	{
		//[ab] - [cd]
		if (a < b)
		{
			if (a >= c && a < d)
			{
				_intr.point = a; _intr.res_type = POINT; return true; 
			}
			if (c >= a && c < b)
			{
				_intr.point = c; _intr.res_type = POINT; return true; 
			}
		}
		else // [ba] - [dc]
		{
			if (a > d && a <=c)
			{
				_intr.point = a; _intr.res_type = POINT; return true; 
			}
			if (c > b && c <=a)
			{
				_intr.point = c; _intr.res_type = POINT; return true; 
			}
		}
	}
	return false; 
}

/**
 * Compute exact intersection of two 2D half-segments. 
 - Half-segment is an oriented segment of which the left side is closed and the right side is open. 
 - In one dimension its form is [x, y). 
 */
template<typename Value>
bool intersect_hseg_hseg_2d (const Segment<Value> & seg1, const Segment<Value> & seg2, IntersectResult<Value> & _intr)
{
	Value  s, t;       /* The two parameters of the parametric eqns. */
	Value num, denom;  /* Numerator and denoninator of equations. */

	const Vector<Value> & a = seg1.p0; 
	const Vector<Value> & b = seg1.p1; 

	const Vector<Value> & c = seg2.p0; 
	const Vector<Value> & d = seg2.p1;

	bool retValue = false;
 
	const int X = 0; 
	const int Y = 1; 
	denom = a[X] * ( d[Y] - c[Y] ) +
		b[X] * ( c[Y] - d[Y] ) +
		d[X] * ( b[Y] - a[Y] ) +
		c[X] * ( a[Y] - b[Y] );

	/* If denom is zero, then segments are parallel: handle separately. */
	if (denom == Value(0.0))
	{
		return intersect_hseg_hseg_2d_parallel(seg1, seg2, _intr);
	}

	num =    a[X] * ( d[Y] - c[Y] ) +
		c[X] * ( a[Y] - d[Y] ) +
		d[X] * ( c[Y] - a[Y] );
	//if ( (num == Value(0.0)) )  retValue = true; 
	//if ( (num == Value(0.0)) || (num == denom) )  retValue = true; 
	s = num / denom;


	num = -( a[X] * ( c[Y] - b[Y] ) +
			 b[X] * ( a[Y] - c[Y] ) +
			 c[X] * ( b[Y] - a[Y] ) );
	//if ( (num == Value(0.0)) )  retValue = true; 
	//if ( (num == Value(0.0)) || (num == denom) )  retValue = true; 
	t = num / denom;
//	if ( (s == 1.0 && t != 0.0) || (s != 0.0  && t == 1.0) ) retValue = true;
	if ( (s == 0.0 && t != 1.0) || (s != 1.0  && t == 0.0) ) retValue = true;

	if      ( (Value(0.0) < s) && (s < Value(1.0)) &&
			  (Value(0.0) < t) && (t < Value(1.0)) )
	{

		retValue = true; 
	}
	else if ( (Value(0.0) > s) || (s > Value(1.0)) ||
			  (Value(0.0) > t) || (t > Value(1.0)) )
	{
		retValue = false;
	}

	_intr.res_type = POINT ;
	_intr.point = Vector<Value> (a[X] + s * ( b[X] - a[X] ), a[Y] + s * ( b[Y] - a[Y] ));

	return retValue;
}


template<typename Value>
bool intersect_convex_poly_poly_2d(const Polygon <Value> & _pol_p, const Polygon <Value> & _pol_q, IntersectResult<Value> & _intr)
{
	unsigned idx_p0 = 0, idx_q0 = 0;
	int count_p = 0, count_q = 0; 
	int size_p = int(_pol_p.size()); 
	int size_q = int(_pol_q.size());
	Vector<Value> origin (Value(0.0), Value(0.0)); /* (0,0) */
	Vector<Value> p0, p1, q0, q1; 
	const int Pin = 1, Qin = 2, Unknown = 0;  /* {Pin, Qin, Unknown}: which inside */
	int inflag = Unknown ;
	int idx_p1, idx_q1; 
	std::vector<Vector<Value> > result ;

	bool first_point = true; 
	int off_p = 0, off_q = 0; 

	IntersectResult<Value> intrseg; 
	bool res = false; 
	Vector<Value> A, B; 
	do 
	{
		idx_p0 = count_p % size_p ;
		idx_q0 = count_q % size_q ; 
		idx_p1 = (idx_p0 + 1) % size_p ;
		idx_q1 = (idx_q0 + 1) % size_q ; 
		A = _pol_p [idx_p1] - _pol_p [idx_p0] ;
		B = _pol_q [idx_q1] - _pol_q [idx_q0] ;

		int cross = orient_2d (origin, A, B);
		int aHB   = orient_2d (_pol_q[idx_q0], _pol_q[idx_q1], _pol_p[idx_p1] );
		int bHA   = orient_2d (_pol_p[idx_p0], _pol_p[idx_p1], _pol_q[idx_q1] );

		std::cout<<"1: "<<Segment<Value>(_pol_p[idx_p0], _pol_p[idx_p1])<<std::endl;
		std::cout<<"2: "<<Segment<Value>(_pol_q[idx_q0], _pol_q[idx_q1])<<std::endl<<std::endl;
		std::cout<<" cross aHB bHA "<<cross<<' '<<aHB <<' '<<bHA<<std::endl;
		Segment<Value > seg1 (_pol_p[idx_p0], _pol_p[idx_p1]);
		Segment<Value > seg2 (_pol_q[idx_q0], _pol_q[idx_q1]); 
		
		res = intersect_hseg_hseg_2d(seg1, seg2, intrseg);
		

		if (res == true )
		{
			if (intrseg.res_type == POINT)
			{
				result.push_back(intrseg.point);
				std::cout<<"Intersect "<<intrseg.point<<std::endl;
				
			} else if (intrseg.res_type == SEGMENT)
			{
				result.push_back(intrseg.segment.p0);
				result.push_back(intrseg.segment.p1);
				std::cout<<"Intersect "<<intrseg.segment.p0<<std::endl;
				std::cout<<"Intersect "<<intrseg.segment.p1<<std::endl;
			}
			if (aHB > 0) inflag = Pin ;
			else if (bHA > 0) inflag = Qin ;
			else inflag = Unknown ;
			
			if (first_point == true)
			{
				first_point = false ;
				off_p = count_p ; //off_p = idx_p1; 
				off_q = count_q ; //off_q = idx_q1;
				std::cout<<"offset "<<idx_p1<<' '<<idx_q1<<std::endl;
			}
		}

		// A & B parallel. Decide its direction. 		
		if (cross == 0 && aHB == 0 && bHA == 0)
		{		
			bool dir_flag = hseg_hseg_2d_parallel_dir(seg1, seg2);
			if (!dir_flag)
			{
				if (inflag == Pin) count_p ++ ;
				else if (inflag == Qin) count_q ++ ;
				else 
				{
					count_p++ ;count_q++;
				}
				inflag = Unknown ;
			}
			else
			{
				if (hseg_hseg_2d_parallel_pqin (seg1,seg2) )
					inflag = Pin; else inflag = Qin;
				if (inflag == Pin) count_p ++ ;
				else if (inflag == Qin) count_q ++ ;
					
			}
		}
		else if (cross >= 0)
		{
			if (bHA > 0)
			{
				if (inflag == Pin && aHB != 0) 
				{
					result.push_back(_pol_p [idx_p1]); 
					std::cout<<"1 Add point "<<_pol_p [idx_p1]<<std::endl; 
				}
				count_p ++ ;
			}
			else 
			{
				if (inflag == Qin && aHB != 0) 
				{
					result.push_back(_pol_q [idx_q1]); 
					std::cout<<"2 Add point "<<_pol_q [idx_q1]<<std::endl;;
				}
				count_q ++ ;
			}
		}
		else //if (cross < 0)
		{
			if (aHB > 0)
			{
				if (inflag == Qin && bHA != 0) 
				{
					result.push_back(_pol_q [idx_q1]);
					std::cout<<"3 Add point "<<_pol_q [idx_q1]<<std::endl; 
				}
				count_q ++ ;
			}
			else 
			{
				if (inflag == Pin && bHA != 0) 
				{
					result.push_back(_pol_p [idx_p1]); 
					std::cout<<"4 Add point "<<_pol_p [idx_p1]<<std::endl; 
				}
				count_p ++ ;
			}
		}
		std::cout<<"count_p "<<count_p<<' ';
		std::cout<<"count_q "<<count_q<<std::endl;
	}
	while ( count_p - count_q < size_p - off_q + off_p && 
			count_q - count_p < size_q - off_p + off_q &&
			(count_p < size_p || count_q < size_q ));

	/* Deal with special cases.  */
	if (result.size() == 0)
	{
		if (is_in_polgyon_2d(_pol_q, _pol_p[0]))
		{
			_intr.polygon  = _pol_p ;
			_intr.res_type = POLYGON; 
			return true; 
		}
		else if (is_in_polgyon_2d(_pol_p, _pol_q[0]))
		{
			_intr.polygon = _pol_q ;
			_intr.res_type = POLYGON ;
			return true; 
		}
		else return false; 
	}

	_intr.res_type = POLYGON ;
	_intr.polygon = Polygon<Value> (result);
	return true ;

}
#endif
